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Simulating open quantum dynamics with time-dependent variational matrix product 
states: Towards microscopic correlation of environment dynamics and reduced system 

evolution 

Florian A. Y. N. SchrodeiQ and Alex W. Chin 
Cavendish Laboratory, University of Cambridge, J. J. Thomson Avenue, Cambridge, CB3 OHE, UK 

We report the development of an efficient many-body algorithm for simulating open quantum 
system dynamics that utilizes a time-dependent variational principle for matrix product states to 
evolve large system-environment states. Capturing all system-environment correlations, we repro¬ 
duce the non-perturbative, quantum-critical dynamics of the zero temperature spin-boson model, 
and then exploit the many-body information to visualize the complete time-frequency spectrum of 
the environmental excitations. Our ’environmental spectra’ reveal correlated vibrational motion 
in polaronic modes which preserve their vibrational coherence during incoherent spin relaxation, 
demonstrating how environment information could yield valuable insights into complex quantum 
dissipative processes. 


I. INTRODUCTION 

Dissipative quantum dynamics can now be probed in 
microscopic, real-time detail, yielding unprecedented in¬ 
sight into system-environment processes whose under¬ 
standing and control will be essential for future quantum 
technologieJi^. Recently, experimental obser vations of 
coherence in organic and biological materialJ^^^IIEHIo!^ 
have strongly motivated a better understanding of the 
microscopic origin and role of ’ultrafast’ (<ps) effects re¬ 
sulting from quantum correlations, memory, bath struc¬ 
ture and non-perturbative system-bath couplings. While 
advanced reduced density matrix techniques can account 
for these phenomengPUOSl^ information is inevitably dis¬ 
carded when tracing out the environment, and a truly 
many-body approach is required for deeper insight into 
the mechanisms at play. However, this necessitates the 
evolution of a macroscopically large system-environment 
state, and the determination of open-system ground 
states and dynamics are only tractable with power¬ 
ful computational techniques, such as exact diagonal- 
ization, multi-configurational Hartree-Fock and various 
time-dependent (density matrix) renormalization group 

technique^i^l*221 

In this article, we present a versatile new approach 
to this problem, based on the recently proposed time- 
dependent variational principle (TDV P) for variational 
matrix product states This many- 

body method gathers together several recent advances 
in VMPS theory to create a fast, efficient algorithm for 
system-bat h dynam ics, where resources can be allocated 
’on the show that the method can cor¬ 

rectly capture the complex, non- Markov ian physics of 
the famous spin-boson model (SBMJPMZl^ and then show 
how visualizing the accompanying environmental dynam¬ 
ics provides an informative, time and frequency-resolved 
spectroscopy of open systems. This combination of ac¬ 
curate system dynamics and powerful diagnostic tools 
for analyzing the detail within system-environment states 
can be applied to a wide range of problems, and could 
be particularly useful for unraveling the physics of the 


’intermediate’ regime of open systemd^. 

The paper is organized as follows. Section |lT] defines 
the model Hamiltonian used in our calculations. Section 
|HI| and |rV| briefly outline the orthogonal polynomial chain 
mapping as well as the VMPS method upon which our 
algorithm is based, including the optimized boson basis 
(OBB) which modifies the matrix product state (MPS) 
network. Section |V| presents the new TDVP scheme re¬ 
quired to time-evolve the modified MPS and presents a 
derivation of the scheme. Finally section [Vl] presents nu¬ 
merical results for the model Hamiltonian and demon¬ 
strates to which detail the environmental dynamics can 
be analyzed and related to system-bath dynamics. 


II. THE SPIN-BOSON MODEL 

The spin-boson model (SBM) has become the bench¬ 
mark for testing advanced open system methods, as well 
as having numer ous d irect applications in physics, chem¬ 
istry and biologji2£122l. R describes a quantum two-level 
system (TLS) that interacts with an environment of har¬ 
monic oscillators via the Hamiltonian [h = 1)P^, 

H = —(Xx— 2 ^ (pri + bli) + 'y tOnblibn, 


where the TLS has an energy bias e, coherent tunneling 
amplitude A and coupling A„ to environmental modes of 
energy uin- The operators ai are Pauli matrices, while 
&]) and bn are bosonic creation and annihilation opera¬ 
tors. Here, we focus on general power-law spectral func¬ 
tions J{uj) = T^Y^n (w ~ = 2TTaiol~‘^UJ^9 (Wc — Oj) 

which parameterize the bandwidth of the environment 
Wc, the (dimensionless) interaction strength a and the 
frequency dependence exponent s which defines sub- 
Ohmic (s < 1), Ohmic (s = 1) and super-Ohmic (s > 1) 
environments. The Ohmic and sub-Ohmic cases possess a 
range of quantum phase transitions ( QPT) ({ (Jz)a.s ^ 0) 
when a exceeds a critical coupling o-JUltlSESHlll 
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Since the VMPS method works especially well on ID 
chain Hamiltonians, it is necessary to transform the 
SBM. We use the orth ogona l polynomial mapping in¬ 
troduced in Prior et aLO^IlSI to obtain a semi-infinite 
ID coupled-chain representation (see next section and 

Ref.li^, 


H = Hs + y CO 


(ao -F aj, j 


L-2 ^ ( 2 ) 

-f [^kalak + tk (4afc+i -h 
k=0 


with effective system-environment coupling cq = 

Jq" ~~ dw and truncation of the chain to a length 

L. 


III. ORTHOGONAL POLYNOMIALS 

The mapping of the environment to a semi-infinite 
chain model is per forme d using orthogonal polynomials 
as described in Ref P^ I '^°i The star-like Hamiltonian 


H = Hs 


g{x)blb^ dx 


^max 

+ A j h[x) (bx + bl) dx, 


( 3 ) 


with the system Hamiltonian Hs and the system opera¬ 
tor A in the interaction term is mapped onto the chain 
Hamiltonian 


H = Hs + Aco 


(ao + aj) 


CXJ 

T ^ ^ -t- tk -t- cik-\-i^k 

k=0 


( 4 ) 


For an arbitrary spectral density J(a;) the mapping can 
be obtained by finding the recurrence relation of polyno¬ 
mials orthogonal with respect to the weight function 


/i2(x) = J(g(x)) 


9'{x) 


( 5 ) 


satisfy the normalization condition 

Xmaa: 

{Pk^Pi)f,= j {x)pk{x)pi{x) dx = 6k^u (7) 
0 

and relate the monic orthogonal polynomials 7rfc(x) via 

7rfc(x) 


Pk{x) = 


Wk{x) 


( 8 ) 


where the norm ||p||^ = is induced by the inner 

product (•, ■)^ under the measure d/i(x) = h^(x)dx. The 
chain parameters co,uJk,tk are related to the coefficients 
Q!fc, Pk of the monic recurrence relation 


TTfc+l (x) = (X - Ofc) TTfc (x) - PkT^k-1 {x) 
LOk = tOcak, 

tk = \/ Pk+lj 

Co = ||7ro{x)||^, 


( 9 ) 


with the initial polynomials 7r_i(x) = 0 and 7ro(x) = 1. 
For the power-law spectral density with hard cut-off at 
the characteristic frequency uJc 


J(io) = 27rQ;Wg (wc — oj ), 


( 10 ) 


the site energies ojk and couplings coAk can be analyti¬ 
cally found as 


I , 

= y (1 


(s + 2k) (2 + s + 2k) J ’ 


— 


(jJc (1 H“ k) (1 -j- 5 H- k) 3 s 2k 


(5 + 2 + 2k) (3 + 5 + 2k) V 1 + 5 + 2A: ’ 


(11) 


for /c = 0 , 1 ,..., L — 2 if the chain is truncated to length 

L. 

This analytic mapping allows a simple inversion = 
J2k Uk{x)a\ to obtain observables of the original Hamil¬ 
tonian from the chain observables. In the presented work 
we used the spin projected displacement 


A/i = , t±(Jzbl + b^ 

Jx \ 2 2 


L-2 


= Hx)Pk{x)^[ d^P'' al)], 


( 12 ) 


k=0 


and the occupation of phonon modes 


where usually the linear dispersion relation g{x) = WcX = 
uj is used without loss of generality. The orthonormal 
polynomials Pk{x) generating the orthogonal transforma¬ 
tion Uk{x) from the continuous variable x onto the dis¬ 
crete chain 

Uk{x) = h{x)pk{x), 

^max 

^k= J Uk{x)bldx, 

0 


L-2 

{bibA = h^{x)pk{x)pi{x) {alai) , (13) 

k,l—0 

where the continuous variable x has to be discretized. 

While for time-evolution the widely used logarithmic 
discretization needs an averaging scheme over multi¬ 
ple calculations to minimize discretization errors, this 
method is exact and gives accurate results in one run. 
In fact this mapping can be r ecover ed in the limit A —>■ 1 
of infinitely fine discretizatiorP^HU. 
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IV. VARIATIONAL MATRIX PRODUCT STATE 
FORMULATION 

We now outline the formalism and features of our al¬ 
gorithm, the intricate details of which can be found in in 

Ref! 17 | 21 | 23 | 

For a ID lattice of size L with sites k and corresponding 
local eigenstate basis \nk) of dimension dk an arbitrary 
state of the Hilbert space \'i') GH can be written as 

dk 

1 ^) = X! |« 1 , • ■ • ( 14 ) 

{rafc} = l 

where the sum is done over every possible combination 
of Ufc. Any |d>) can be written as a matrix product state 
(MPS) I'^MPs) € M C 'H via iterative application of sin¬ 
gular value decompositions (SVD) on resulting 

in the rank-3 tensors A{k) G £Dk-ixDkXdk 

| 4 rMPs)= E ( 15 ) 

{nfe} = l 

where we have open boundary conditions (og = = 

1), use the indexing = {A(k))ak_,,ak,nk = 

{A'^>^{k))ak_^,ak^ and allow to omit the site argument for 
clarity. 

Furthermore we employ an optimized boson basis 
(OBB), as introduced by Guo et al^, which is realized 
via an additional map (isometry) V G C‘^i‘^‘^oBB,k Rom 
the optimized basis jn^) into the local basis \nk) 

doBB,k 

KUak= E ^aUakVnk,i.k. ( 16 ) 

hk = l 

where we will write A instead of A throughout this work. 
This mapping allows high compression {doBB,k ^ dfc) 
of the local oscillator basis in case of large variances 
Yai {uk), which has been shown to be highly effective 
in dealing with quantum critical SBMji3_ 

Once we represented a state as an MPS, any variational 
optimization and time-evolution is performed iteratively 
by sweeping along the chain. During the sweeping proce¬ 
dure we keep the state in a mixed canonical form to ben¬ 
efit from orthogonality conditions. Furthermore only one 
matrix can be focused (centered), which means its vec- 
torization v has unit length {v\v) = 1 = ('I'mpsI'I'mps)- 
Fig. diagrammatically explains all normalizations of A 
and the possible center matrices Ac{k),Vc{k) and bond 
centers Ca{k), Cn{k) which have special relevance for the 
time-evolution explained later. When focusing on site fc, 
all A matrices of site i < k will be kept left-orthonormal 
while matrices i > k will be right-normalized to produce 
orthonormal left and right basis states 14)^), 

|d/(A)) = 

E (17) 

a,k-l,nk,ak 







FIG. 1. (color online), (a)-(c) All possible normalizations of 
A. The straight line represents the identity matrix, (d)-(g) 
Center matrices (orange) and corresponding network normal¬ 
ization for one site. 


The efficiency of MPS is based on low-rank tensor ap¬ 
proximations which significantly reduce the number of 
variational parameters in the tensors A{k),V{k). This 
restricts the Hilbert space spanned by the MPS to a 
manifold Ai G %. The combination of high truncation 
and dynamical ’on the fly’ bond adjustment, while rep¬ 
resenting an optimal manifold allows the implementation 
of an efficient variational algorithm with sign ificant ad- 
vantages for computational speed and accuracjff^^^SEIIIIl 
We truncate and expand the bonds and fik such that 
the smallest singular values kept are within 10“^ and 
10-4.5 leading iq adaptive dimensions Dk and doBB,k 
with upper bounds D^ax and doBB,max- 


V. TIME-DEPENDENT VARIATIONAL 
PRINCIPLE 

The time evolution is performed with the time- 
depe nden t variational principle (TDVP) described in 
HefsjUUll 11 ig derived from the Dirac-Frenkel variational 
principle and obtained by projecting the Schrodinger 
equation onto the tangent space of the MPS manifold 
Ai to find optimal equations of motion for each MPS 
center tensor within Ai to generate the best approxima¬ 
tion l'I'MPs(i)) to the exact state |'I>(t)). Haegeman et 
al. have shown that this is equivalent to a Lie-Trotter 
splitting of the MPS instead of the time-evolution oper- 
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ator U{t) = Unlike the Suzuki-Trotter splitting of 

errors only accrue from the integration scheme. 

Since the OBB modifies the MBS network, our im¬ 
plementation of the TDVP needed an extension of the 
original TDVP scheme. To prove that this can be done 
and to derive the appropriate integration scheme, we 
will present in Sec. |V A| the derivation of the tangent 
space projector PT^^^pg)M in the formalism established 
by Ref.l^. Additionally we will give a more simplified de¬ 


scription of the resulting integration scheme in Sec. VB 


A. Derivation of the tangent space projector 

This derivation will be similar to Theorem 3.1 of ReP^ 
but using the left canonical form of the MPS with OBB 
matrices. To be more consistent with their notation and 
for clarity, we will use n^, hfc and Ofc instead of dk, doBB,k 
and Dk to address bond dimensions. 

Additionally to the definitions of Lubich et al. we need 
to introduce further notation to address the components 
of an MPS tensor X € ^nd their unfoldings 

(as depicted in Fig. [^. 

Any site tensor Ck G jjafc-ixnfcxat written 

as a matrix using the left and right unfolding G 
]j(afc_i,nfc)xafc g ]j(afc,nfc)X^ Similarly the 

orthonormalized site tensors Q^jQk ^ are 

unfolded as 


Q< g ]^(afe-l")s)xafc 

Q> (Z ]^(“fc"-fc)xafc_i 


(18) 


and additionally fulfill the normalization conditions 

Qfe Qk ~ = lafc_i- 

Unfolded segments of the MPS to the left and right 
of site k are denoted as X<k-i G and 

X>k+i G respectively. If these parts are 

left and right orthonormalized we will write Q<k-i and 

This notation can be transferred to the OBB tensors 
Ak G and Vk G forming a decom¬ 

position of the site tensors Ck as 

C< = {Vk®ta,.,)A<. (19) 

The OBB tensors can also be orthonormalized, de¬ 
noted as QA,k and Qv,k G where QYkQv,k = 


A,k — {Qv,k ® Q<k — l^ A.k^>k-\-l * 

= |Ten(fe) i5QKfc(Q^ fc)^(Ai>fe+i 0 Q<k-i)^ 


with X'^> L -I- 1 = 1 since ol = 1. Since these subspaces 


Ql’^QXk = compose and as: 

Qk ~ {Qv,k ® lofc_i) QA,k’ / 2 Q'\ 

Qk ~ {Qv,k ® lafc) QA,k- 

Furthermore it is necessary to define an new unfolding 
G to support the following notation. 

Let be the fcth tensor transpose cycling the fcth 
tensor dimension of X to the 1st position while keeping 
the order of all other dimensions. This permutation can 
be written as 


A(*l U2 j ■ ■ ■ j *l) — A (*cr,, (1) , (2) I • ■ ■ I iak{L))- (21) 

with the fc-cycle Uk = {k 1 2 • • • fc — 1) G S/,. The 
inverse of this transpose is defined by the fc-cycle 

G-k = = {k k — \ ■ ■ ■ 2 1). The fcth unfolding 

X^^'> of X denoted with round braces is then equivalent 
to the 1st unfolding of the fcth tensor transpose of X 


(fe) _ g |^nfcX(ni...nfc_infc+i...ni) 

AT-k (22) 


A = Ten 


(fe) 


x(fc) 


= Teni 


X'^kA) 


with the corresponding tensor reconstruction as its in¬ 
verse. 

A matrix A G > n with rank(A) = n has a 

left-inverse 

A-i = {A^ 

A-^A = In, 

and a projector Pa onto the range of A 

Pa = AA-\ (24) 



which will be used for the definition of the tangent space 
projector. 

The tangent space TxM of any state A G AI of the 
OBB-MPS manifold can be constructed from the orthog¬ 
onal subspaces VA,k,Vv,k- We will use the MPS in left 
canonical gauge including site fc 


A = Teufc 


{Qv,k ® Q<k—l) QA,k^>k +1 


(25) 


Thus we have 


< , G = 0 for fc ^ l} , 

■ ■.6Q<,eR^'‘^^\QXk5Qv,k=0}, 


(26) 


are mutually disjoint, the tangent space is 

L 

TxM. = {^^A,k ® Tv,k), 

fc=i 


( 27 ) 

















5 


and allows the decomposition of any tangent vector 6X G 
TxM 


5X = ^ {5XA,k + SXv,k ), , , 

fc=i 1 

(SXi^kySXj^i) =0 i0Ti^j,k^l. 


where we used the projectors 

Pv,k = Qv,kQv,k^ 

P<k — 1 — Q<k—lQ<k — l^ 

P>k+1 = -^>fc+l-^>fe+l “ Q>fc+lQ>fc+l- 


(36) 


We want to find the tangent space projector PtxM 

for an arbitrary Z G - such that the projec- In the notation of Reff^ the tangent space projector 

tion SU = Ptxm(^) S PxM. is orthogonal and therefore thus reads: 

satisfies 


{6U,dX) = {Z,SX) y6X G TxM. (29) 
Due to the disjoint subspaces Vi^k we can decompose 


L 

6U = J2iSUA,k + SUv,k), (30) 

with SUi^k G Vi,fc and find 

{SU,,k,SX,,k) = {Z,6X,^k), (31) 

for i G {A,V}- The projection is then defined by the SB 
matrices in 


Sui% = iQv,k 0 Q<k-i)SB<,X^,^,, 
5U\^1 = ,)^(X>fe+i (8) Q<k-if, 


with Q^\5B^ ^=Q ioT k ^ L and Qyj^SBy.k = 0. 

To find expressions for 5B“^ f. and SBv,k we follow the 
same steps as RefP. Here we give an outline for 6Bv,k to 


obtain SUyl. First substitute the expressions for SUyj. 

and isolate the term 


(fe) 


and SXl^l (Eq. 


26) into Eq. 


SQv,k in both inner products. 
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{sul,%sxl,^l) = {SBv,kM^,SQv.kM^) 

= {dBv^kM^M,SQv,k), (33) 
(ZW,JX|.5) = {Z('^'>M,6Qv,k), 

where we defined M = (X>fc+i(8)Q</i:-i)Q^ f. for clarity. 
To remove the inner products on both sides we use the 
property {A, C) = (B, C) -G Pc A = PqB with ^ = 
In;, - Pv,k and then substitute the found expression for 
SBv,k into SUH^I 


SBv.k = (In, - Pv,fc)Z«M(M^M)-\ 
SUH^l = (1„, - Pv,k) Z('=)(MM-i), 


Thus we derived the projection of Z onto the tangent 
space as 


6U[% = [{Pv,k ® P<k-i) - P<k] Z^^'^P>k+i. 

= (1„, - Py.fc) zWPjvf. 


MPs) Z-^ ^ L 




= E^J 

k^l 

k^l 

L 


(37) 


(Ife-P?')- 


k^l 


with the left and right chain projectors 


p[l:fc] 


p[k+l\L] 

Ft 


P 


[fc] 


A[l:fc-l,fc+l:L] 

P4 


t h-K!>( 4 -K!i, 

a, = l 

ak — l 
dOBB,k 

T. i*S>(fgi. 

rik — l 
doBB,k 

nfc = l 

(38) 


Using these expressions to project H\'i>MPs) onto 
Ti^mps)-^ we get, after a Lie-Trotter decomposition of 
the MRS, a system of equations of motion for each MRS 
tensor Vc{k) and Ac{k) as well as for the bond matrices 
Cn{k) and Ca{k) which can be solved directly by 


Fc(fc,i) = e-*^'^('=)Wc(fc,0), 
Cf,(fc,t) = e+‘""W‘Cf;(fc,0), 
Ac(fc,t) = e-*"-^('=)‘Ac(fc,0), 
C„(fc,t) = e+*"“W‘CJA:,0). 


The effective local Hamiltonians of each tensor used in 
the exponentials are obtained by a full contraction of the 
Hamiltonian with the MRS, while omitting the respective 
tensor. A diagrammatic representation is given in Fig. 
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(c) Hv(k) 



(e) Hf,ik) 


FIG. 2. All effective local Hamiltonians needed for the TDVP 
with OBB. 


B. l-tensor update scheme 

A key formal result we present is this new OBB-TDVP 
scheme with a single-tensor update, which is given by 

Vcik, t + At) = t), 

iSVD 

i SVD (39) 

Ac{k,t + At) = 

iSVD 

Ca{k,t) = t + At), 


where bold notation indicates vectorizations of the ten¬ 
sors. The effective local Hamiltonian H^{k) generating 
the optimal time-evolution of a center matrix is obtained 
by full contraction with all MBS matrices except for the 
to be evolved center matrix as depicted in Fig. This 
scheme can be seen to maintain the general form of the 
TDVP equations given in Refs.^^HlI! Essentially, after 
each forward evolution t —> t-|-At of a centered tensor Vc 
or Ac, the centered bond matrix C^jJ^t -|- At) obtained 
via SVD has to be evolved backwards in time t -I- At —>■ t 
before contraction with the next tensor. By sweeping 
left-to-right and back to the left with a timestep At/2 
one obtains a second order symmetric integrator with er¬ 
ror of order 0{At^). Since for each evolution step the 
entire effective Hamiltonian is applied, it is possible to 
include long-range interactions. 

The single-tensor update with OBB has signifi¬ 
cant advantages over the 1-site update TDVP with¬ 
out OBB in cases where a large number of local 
states has to be considered. The two additional evo¬ 
lution steps for Vc and Cn can outweigh the un¬ 
favorable scaling of the 1-site TDVP without OBB. 
The computational complexity of the TDVP with OBB 
scales as 0{max{D^doBB, D^dQBB>doBBd^,dQg^d)) 
which is significantly faster than the original scaling 
0{max{D^d, D'^d^)) without OBB for doBB ^ d. The 
scaling behavior of OBB-TDVP is also more prefer¬ 
able in comparison with the common time-evolving 
block decimation (TEBD) {0{d^D^)) and the very re¬ 
cently proposed TEBD with local basis optimization 
(TEBD-LBO) of Brockt et al^ with a scaling of 
0{u\ax{d^BB^^ t^^D"^)). The latter uses a technique 
very similar to the OBB to map the local Hilbert space 
onto an optimized basis to reduce the computational ef¬ 
fort. 

As evident from the example given in Fig. |^it is pos¬ 
sible to increase the local Hilbert space on each site from 
dk = 200 to dk = 500 at no extra cost by introducing 
an OBB of doBB = 65. Similarly a calculation need¬ 
ing dk = 500 can be accelerated by a factor of 4-10 by 
using an OBB, depending on the amount of (allowed) 
entanglement. Generally, the smaller doBB, the less en¬ 
tanglement can be captured by the MPS between a site 
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FIG. 3. (color online). Comparison of CPU time per sweep 
for TDVP without OBB and = 200, dk = 500 (lines) and 
for TDVP with OBB and dk = 500 (crosses). 

and its left and right environment, but the faster and 
more memory efficient is the computation. 

VI. SIMULATED SBM DYNAMICS 

If not otherwise stated, we take A = 0.1, e = 0,a;c = 
1, At = 0.1 — 1 with initial MPS dimensions dk = 30, D = 
5,doBB = 5 and maximal Dmax = ^,doBB,max = 15 
in all simulations. All results had converged sufficiently 
w.r.t. U by H = 5, as shown in the appendix. The 
required chain length L depends on Wc, the simulated 
time range T and the extracted observables. If only the 
system dynamics are desired, we take L = '^T, while for 
environment observables in frequency space, we need L = 
to avoid artifacts caused by unphysical reflections at 
the end of the chain (see Fig. |^. For reference, a single 
sweep with dimensions dk = 30,1? = 5,doBB = 15, and 
L = 200 takes 4 seconds on one core of an Intel Core 
17-4790 CPU. 

We will consider two different initial state preparations 
for time-evolution, both with the environment at zero 
temperature. The polarized coupled state is obtained via 
variational ground state optimization of the z-pol arized 
spin coupled to the environment, as in Refs.ll^H^. The 
uncorrelated product state has the vacuum state |f2) of 
the bath coupled to the z-polarized spin. 

A. Spin-dynamics 

In Figs. |4]|^ we present the dynamics in the weak 
{a < ac) and strong coupling regime {a > ac) for Ohmic 
s = 1 and sub-Ohmic s = 0.25,0.5,0.75 spectral densi¬ 
ties. Globally, we find excellent agreement with previous 
numerical work on the non-Markovian dynamics of the 
SBM, such as the time-dependent numerical renormal¬ 
ization group (TD-NRG)(13and path integral methodJi^. 
In brief, the Ohmic results (Fig. show increasingly 


0 50 100 150 200 250 300 



FIG. 4. (color online), (a) Weak {a < Oc) and (b) strong cou¬ 
pling for Ohmic (s = 1, Oc = 1) spectral density at different 
a. (b) was obtained at At = 0.01. 



FIG. 5. (color online), (a)-l-(b) Weak and (c)-l-(d) strong 
coupling (below and above a^) for sub-Ohmic spectral density 
(s = 0.5, Oc = 0.107) at different couplings a. Initial state is 
(a)-l-(c) the product state and (b) + (d) the coupled state. 


damped oscillations as a increases for a < 0.5, over¬ 
damped relaxation for 0.5 < a < 1 which relax more 
slowly as a —t 1 and complete localization above the 
quantum critical coupling of ttc = 1- 

The spin dynamics of the sub-Ohmic SBM close to 
the coherent-incoherent changeover aci{s) are shown 
under the two different initial preparations of product 
state and coupled state in Fig. and In con¬ 
trast to the Ohmic case, the sub-Ohmic dynamics al¬ 
ways remain underdamped for s = 0.25, showing at 
least one osc illation even above = 0.022, as re¬ 
cently found Furthermore an initial polarization 

of the bath leads to a higher frequency of spin oscilla¬ 
tions with stronger coupling a which is not observable 
in the product state since the polarization persists on 
a larger time scale, giving an effective bias to the TLS 
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LUct UJct 

FIG. 6. (color online). Weak and strong coupling (below 
and above a^) for sub-Ohmic spectral densities ((a) + (b): s = 
0.25, ac = 0.022; (c) + (d): s = 0.75, « 0.3, aci ^ 0.22) at 

different couplings a. Initial state is the product (a) + (c), or 
the coupled state (b) + (d). 


HamiltoniarP^.The frequency of these oscillations is non¬ 
monotonic, initially decreasing due to dressing of the tun¬ 
neling matrix element (a system-bath correlation effect) 
and then increasing for stronger coupling. The s = 0.5 
case exhibits overdamping only with the product state for 
a > aci ~ etc = 0.107, while the coupled state always 
has initial oscillations even at strong coupling. The final 
value of (cTj) reflects the mean field nature of the QPT 
in systems with 0 < s < 0.5, for which the spin magne¬ 
tization growscontinuously from zero above the critical 
couplin^^^i^^^ISIl g _ Q 75 both the product state 
and the coupled state lead to overdamped dynamics for 
a ^ aci ~ 0.22 and localize for a > ac ~ 0.3. The ac¬ 
curacy of the calculations across a wide range of spectral 
densities from weak to strong coupling combined with the 
efficient use of computational resources demonstrates the 
versatility of the VMPS implementation. 


B. Environment dynamics and spectroscopy 


Having verified our method, we now use an efficient 
inversion of the chain mapping (see Sec. Ill) to present 
the dynamics of the entire environment in time-frequency 
spac^^. Fig. [^a) and its inset shows the population 
of each mode for intermediate Ohmic coupling, with an 
initially broad excitation and subsequent emergence of a 
sharp resonance peak around the TLS energy gap. The 
peak in this environmental ’absorption spectrum’ rises 
on the timescale of spin relaxation (uct « 300), while its 
position evolves and is complete by ujet « 110. The final 
peak position is reached after about one period of the 
resonant frequencies of the environment modes {uires ~ 
O.OOwc), consistent with the ’sampling’ required for the 
broadband environment to resolve the TLS gap. A range 


of novel phenomena related to non-detailed balance and 
phase-dependent relaxation have recently been predicted 
to exist prior to this tim^^, though we will not explore 
this further here. 

Instead, we note that this timescale maybe addition¬ 
ally modified by polaronic dressing (spin-bath entangle¬ 
ment), which leads to a renormalization (suppression) 
of the bare TLS energy gap A that drives Ohmic and 
sub-Ohmic ground states towards their QPTs. The envi¬ 
ronment spectrum in Fig. j^a) clearly resolves the emer¬ 
gence of the renormalized TLS’s energy gap A^. This 
ultrafast process, dominated by high frequency (wk > A) 
modes, is generally hard to observe but important in 
organic exciton transfer and has been seen in inorganic 
semiconductor^. According to Silbey and Harris’ varia¬ 
tional polaron theory for the ground state of the Ohmic 

SBAPI, A,. = A (^ which agrees with the peak 
position extracted from the environmental spectra. This 
is shown in Fig. j^b) for s = 1 and s = 3 in the regime 
0 < a < 1. In the chain basis the renormalization is ac¬ 
companied by the persistent excitation of the sites closest 
to the system, which defines an effectively screened sys¬ 
tem ’seen’ by the rest of the environment (see Fig. 

- an observation familiar in NRG studies of the related 
Hondo problenP^. The dynamics of the collective coor¬ 
dinate of the first chain site were also used to verify a 
novel coherence pumping mechanism in a recent study of 
a simple photosynthetic pigment-protein complejP^. 

Accompanying the emission of phonons into the res¬ 
onant modes, we also observe prominent, damped oscil¬ 
lations for all modes, which vanish on the timescale of 
spin relaxation (Fig. [^. Recently, Hera et al. have 
shown that strong coupling induces intermode entangle¬ 
ment which could i nduce apparent mixing and damping 
of individual mode^^HM!^ this is predicted for slow 
and resonant modes, whereas we find damped oscillations 
at all frequencies. Given the unitary, energy preserving 
evolution of the many-body state, this observation war¬ 
rants further consideration. To investigate this, we per¬ 
form a type of state-selective coherent spectroscopy to 
look at further details within the many-body state. Fig. 
shows the displacement of modes projected onto 

the up spin state, i.e. we compute ). 

It is immediately seen that the displacements do not show 
any damping and oscillate at their natural frequencies 
over the length of the simulation. 

Figure shows the typical dynamics of a high fre¬ 
quency mode. The positive and negative displacements 
for I t) and | ),) are characteristic of polaronic entangle¬ 
ment between the spin and oscillator. However, it can 
be seen that their motion is highly correlated; maximum 
mode displacement on one spin state is always obtained 
at the minimum of the other. Moreover, the momen¬ 
tum of the mode in each spin state is the same and pre¬ 
served; it is not randomized by the dissipative spin-flip 
dynamics. This classical correlation within the entangled 
state is important, and explains both the appearance of a 
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FIG. 7. (color online), (a) Phonon modes uJk occupation for an Ohmic environment with a = 0.2 exhibit a strong resonance 
at the renormalized tunneling amplitude Ar « 0.057 and a damping of oscillations, (b) Simulated renormalized tunneling Ar 
versus analytic expression (lines) shows good agreement. Polaron theory predicts A,- = Ae~“^ for s = 3. (c) The projected 
phonon displacement fl for s = 1, a = 0.4 oscillates even after spin relaxation at t « 350. 


logio ("l) 



20 40 60 80 100 120 140 160 180 200 

Site k 


FIG. 8. (color online). The occupation of the chain for 
s = 1 , 0 ? = 0.2 shows a clear constant offset on the first few 
sites indicating the dressing and renormalization of the spin. 
Additionally the chain exhibits reflections from the end at 
uict ~ 380 leading to artifacts in the inverse chain mapping. 
The open system dynamics are correct until the arrival of the 
reflected waves at ~ 760. 


stationary renormalized and the apparent relaxation 
of the mode populations {b\hk)- To motivate this, con¬ 
sider the trial wavefunction |'I'(t)) = C^{t) |t) \flit)) + 
li) l/fe (0)? where we neglect the other environ¬ 
mental modes. The oscillator wavefunctions for each 
spin state are time-dependent coherent states \fl{t)) = 
e^fkW^l-fk it)bk) |Q^_ expectation value (ax) is then 
(ax) = 2'Si[C^(t)Cl{t) {ftit)\fl(t))]. The oscillatory dis¬ 
placements we find fit 2fl{t) « (1 — e^*"*) and 

2f^{t) « gfeW^^(l-|-e*“''‘^), after about one oscillation pe- 



0 100 200 300 400 500 


UJct 

FIG. 9. (color online), (a) The occupation (rifc) and (b) the 
spin-projected displacements of the mode oJk = 0.325 
with initial product state. 



Mode ujk /iOc 

FIG. 10. (color online). The phonon occupation for s = 
0.75, a = 0.4 with initial product state oscillates even after 
the relaxation of the spin oscillations at uj^t « 150 due to the 
finite (ctz) yf 0. 
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riod. As {ax) is determined by the overlap of the oscilla¬ 
tor wavefunctions entangled with each spin state, we see 
that the relative displacement of these wavefunctions is 

nL "T* 1 

preserved for all further times, obeying f{^ — fl = gk^^k ■ 
Thus the correlated motion provides a constant renor¬ 
malization of the spin tunneling, which when summed 
over all modes leads to the predicted by ground state 
theories and observed here as the peak position of the 
environmental absorption. Indeed, had the relative dis¬ 
placement been time-dependent then {ax) would not re¬ 
lax to a stationary value. Next we see that the mode 
population {hlhk) = \C^{t)\'^\fl{t)? + |C';WPl/fc(t)P = 
\gk^k^{^ -I- {a^) cos{uJkt)). Again, we see the correlated 
motion only leads to oscillations of the population when 
the spin is out of equilibrium {{ax) 1), and explains 

why the apparent population damping for all high fre¬ 
quency modes is set by the relaxation of the spin. In¬ 
terestingly, this analysis predicts persistent population 
oscillations for a > Oc in sub-Ohmic baths, or a biased 
TLS, which is confirmed in Fig. |10[ However, we note 
that in real systems the environment must be connected, 
albeit weakly, to an external bath, and oscillations will 
vanish at very long times {all observables become sta¬ 
tionary) . 

These results show how broadband environmental dy¬ 
namics can be analyzed at their various significant fre¬ 
quency scales. Indeed, we have only considered a small 
subset of the SBM physics; similar analysis must also be 


applied to resonant and slow frequency modes, the dy¬ 
namical hierarchy and feedback between their contribu¬ 
tions and the role of any isolated modes, to build the full 
picture. This dynamical dissection, made possible by our 
VMPS many-body approach, not only offers new insights 
into the fundamental aspects of complex open systems, 
it may also guide the construction of cheape r ansdtze for 
describing their ground states and dynamicJii^^^IlSIllIlI! 
More immediately, the observation of conserved vibra¬ 
tional coherences during heavily damped spin relaxation 
has relevance for time-resolved observation of transfer, 
activation and deactivation of vibrational coherences, re¬ 
cently shown to be a powerful tool for exploring ultrafast, 
coherent processe s in optoe lectronic systems and biologi¬ 
cal photoreaction^SEfillHIlll, Finally, we note that further 
extensions are required for general applications, many of 
which - such as finite temperatures and damping of the 
primary environment - have already been developed in 
related methods, such as the DMRG-based TEDOPA al¬ 
gorithm of Prior et al^. 
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